This script analyzes and plots data for Symbiont Integration 2020 respirometry data. Plots are displayed in Plotting section. Results are provided in Analysis section along with summaries of potential implications.
Set up workspace, set options, and load required packages.
knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE)
## install packages if you dont already have them in your library
if ("tidyverse" %in% rownames(installed.packages()) == 'FALSE') install.packages('tidyverse')
if ("car" %in% rownames(installed.packages()) == 'FALSE') install.packages('car')
if ("lme4" %in% rownames(installed.packages()) == 'FALSE') install.packages('lme4')
if ("lmerTest" %in% rownames(installed.packages()) == 'FALSE') install.packages('lmerTest')
if ("scales" %in% rownames(installed.packages()) == 'FALSE') install.packages('scales')
if ("cowplot" %in% rownames(installed.packages()) == 'FALSE') install.packages('cowplot')
if ("ggplot2" %in% rownames(installed.packages()) == 'FALSE') install.packages('ggplot2')
if ("effects" %in% rownames(installed.packages()) == 'FALSE') install.packages('effects')
if ("emmeans" %in% rownames(installed.packages()) == 'FALSE') install.packages('emmeans')
if ("multcomp" %in% rownames(installed.packages()) == 'FALSE') install.packages('multcomp')
#load packages
library("ggplot2")
library("tidyverse")
library('car')
library('lme4')
library('lmerTest')
library('scales')
library('cowplot')
library('effects')
library('emmeans')
library('multcomp')
Load data from LoLinR.
PRdata<-read.csv("Mcap2020/Output/Respiration/oxygen_P_R_calc.csv") #load data
Separate project specific data.
#remove all rows of wells that did not have samples or blanks
PRdata<-PRdata[!is.na(PRdata$Type),]
#format columns
PRdata$dpf<-as.factor(PRdata$dpf)
#subset fused juvenile data frame for a different project, then save as separate file
fused_oxygen <- PRdata[which(PRdata$Lifestage=='Juvenile'),]
fused_oxygen <- droplevels(fused_oxygen)
Calculate a P:R ratio using gross photosynthesis.
PRdata$ratio<-abs(PRdata$GP.nmol.org.min)/abs(PRdata$R.nmol.org.min) #calculate ratio with absolute values
#remove outliers detected by values of P:R ratio data
PRdata<-PRdata%>%filter(ratio < 15)
fused_oxygen$ratio<-abs(fused_oxygen$GP.nmol.org.min)/abs(fused_oxygen$R.nmol.org.min)
Load sample metadata to identify stage and hpf.
metadata<-read.csv("Mcap2020/Data/Respiration/PR_metadata.csv")
PRdf<-PRdata%>%
filter(!Lifestage=="Juvenile")
#load in sample information
PRdf$code<-paste(PRdf$Date, PRdf$Lifestage)
metadata$code<-paste(metadata$Date, metadata$Lifestage)
PRdf$hpf<-metadata$hpf[match(PRdf$code, metadata$code)]
PRdf$group<-metadata$group[match(PRdf$code, metadata$code)]
PRdf$hpf<-as.factor(PRdf$hpf)
Plot data as a scatterplot
r_plot<-PRdf %>%
ggplot(., aes(x = hpf, y = R.nmol.org.min)) +
geom_hline(yintercept=0, linetype="dashed", color="black", size=0.75)+
geom_boxplot(aes(color=Lifestage, group=group), outlier.size = 0, position = position_dodge(1), lwd=1) +
#geom_smooth(method="loess", se=TRUE, fullrange=TRUE, level=0.95, color="black") +
geom_point(aes(fill=Lifestage, group=group), pch = 21, size=4, position = position_jitterdodge(0.2)) +
xlab("Hours Post-Fertilization") +
scale_fill_manual(name="Lifestage", values=c("#DFC27D", "#80CDC1","#003C30"))+
scale_color_manual(name="Lifestage", values=c("#DFC27D", "#80CDC1","#003C30"), labels=c("Embryo", "Larvae", "Metamorphosed \nRecruit"))+
ylab(expression(bold(paste("R (nmol ", O[2], " larva"^-1, "min"^-1, ")")))) +
scale_y_continuous(limits=c(-0.1, 0.1), labels = scales::number_format(accuracy = 0.01, decimal.mark = '.'))+
theme_classic() +
theme(
legend.position="none",
axis.title=element_text(face="bold", size=16),
axis.text=element_text(size=12, color="black"),
legend.title=element_text(face="bold", size=14),
legend.text=element_text(size=12)
); r_plot
#EGG: #8C510A
#EMBRYO: #DFC27D
#LARVAE: #80CDC1
#RECRUIT: #003C30
#ggsave("Mcap2020/Figures/Respiration/Respiration.png", r_plot, dpi=300, w=8.5, h=5, units="in")
Plot data as a scatterplot
p_plot<-PRdf %>%
ggplot(., aes(x = hpf, y = P.nmol.org.min)) +
geom_hline(yintercept=0, linetype="dashed", color="black", size=0.75)+
geom_boxplot(aes(color=Lifestage, group=group), outlier.size = 0, position = position_dodge(1), lwd=1) +
#geom_smooth(method="loess", se=TRUE, fullrange=TRUE, level=0.95, color="black") +
geom_point(aes(fill=Lifestage, group=group), pch = 21, size=4, position = position_jitterdodge(0.2)) +
xlab("Hours Post-Fertilization") +
scale_fill_manual(name="Lifestage", values=c("#DFC27D", "#80CDC1","#003C30"))+
scale_color_manual(name="Lifestage", values=c("#DFC27D", "#80CDC1","#003C30"), labels=c("Embryo", "Larvae", "Metamorphosed \nRecruit"))+
ylab(expression(bold(paste("P (nmol ", O[2], " larva"^-1, "min"^-1, ")")))) +
scale_y_continuous(limits=c(-0.1, 0.1), labels = scales::number_format(accuracy = 0.01, decimal.mark = '.'))+
theme_classic() +
theme(
legend.position="none",
axis.title=element_text(face="bold", size=16),
axis.text=element_text(size=12, color="black"),
legend.title=element_text(face="bold", size=14),
legend.text=element_text(size=12)
); p_plot
#EGG: #8C510A
#EMBRYO: #DFC27D
#LARVAE: #80CDC1
#RECRUIT: #003C30
#ggsave("Mcap2020/Figures/Respiration/NetPhotosynthesis.png", rp_plot, dpi=300, w=8.5, h=5, units="in")
Plot data as a scatterplot
gp_plot<-PRdf %>%
ggplot(., aes(x = hpf, y = GP.nmol.org.min)) +
geom_hline(yintercept=0, linetype="dashed", color="black", size=0.75)+
geom_boxplot(aes(color=Lifestage, group=group), outlier.size = 0, position = position_dodge(1), lwd=1) +
#geom_smooth(method="loess", se=TRUE, fullrange=TRUE, level=0.95, color="black") +
geom_point(aes(fill=Lifestage, group=group), pch = 21, size=4, position=position_jitterdodge(0.2)) +
xlab("Hours Post-Fertilization") +
scale_fill_manual(name="Lifestage", values=c("#DFC27D", "#80CDC1","#003C30"))+
scale_color_manual(name="Lifestage", values=c("#DFC27D", "#80CDC1","#003C30"), labels=c("Embryo", "Larvae", "Metamorphosed \nRecruit"))+
ylab(expression(bold(paste("GP (nmol ", O[2], " larva"^-1, "min"^-1, ")")))) +
scale_y_continuous(limits=c(-0.1, 0.1), labels = scales::number_format(accuracy = 0.01, decimal.mark = '.'))+
theme_classic() +
theme(
legend.position="none",
axis.title=element_text(face="bold", size=16),
axis.text=element_text(size=12, color="black"),
legend.title=element_text(face="bold", size=14),
legend.text=element_text(size=12)
); gp_plot
#EGG: #8C510A
#EMBRYO: #DFC27D
#LARVAE: #80CDC1
#RECRUIT: #003C30
#ggsave("Mcap2020/Figures/Respiration/NetPhotosynthesis.png", gp_plot, dpi=300, w=8.5, h=5, units="in")
Plot data as a scatterplot
pr_plot<-PRdf %>%
ggplot(., aes(x = hpf, y = ratio)) +
geom_hline(yintercept=1, linetype="dashed", color="black", size=0.75)+
geom_boxplot(aes(color=Lifestage, group=group), outlier.size = 0, position = position_dodge(1), lwd=1) +
#geom_smooth(method="loess", se=TRUE, fullrange=TRUE, level=0.95, color="black") +
geom_point(aes(fill=Lifestage, group=group), pch = 21, size=4, position = position_jitterdodge(0.2)) +
xlab("Hours Post-Fertilization") +
scale_fill_manual(name="Lifestage", values=c("#DFC27D", "#80CDC1","#003C30"), labels=c("Embryo", "Larvae", "Metamorphosed \nRecruit"))+
scale_color_manual(name="Lifestage", values=c("#DFC27D", "#80CDC1","#003C30"), labels=c("Embryo", "Larvae", "Metamorphosed \nRecruit"))+
ylab(expression(bold(paste("P:R")))) +
scale_y_continuous(limits=c(-0.5, 10), labels = scales::number_format(accuracy = 0.01, decimal.mark = '.'))+
theme_classic() +
theme(
legend.position="right",
axis.title=element_text(face="bold", size=16),
axis.text=element_text(size=12, color="black"),
legend.title=element_text(face="bold", size=16),
legend.text=element_text(size=14)
); pr_plot
#EGG: #8C510A
#EMBRYO: #DFC27D
#LARVAE: #80CDC1
#RECRUIT: #003C30
#ggsave("Mcap2020/Figures/Respiration/NetPhotosynthesis.png", gp_plot, dpi=300, w=8.5, h=5, units="in")
full_fig<-plot_grid(r_plot, p_plot, gp_plot, pr_plot, ncol=4, nrow=1, rel_heights= c(1,1,1,1), rel_widths = c(1,1,1,1.2), label_y=1, align="h")
ggsave(filename="Mcap2020/Figures/Respiration/respirometry_fig.png", plot=full_fig, dpi=500, width=26, height=6, units="in")
Build linear mixed effect model and examine for Respiration.
Rmodel1<-lmer(R.nmol.org.min~group + (1|Run), data=PRdf) #run as random
anova(Rmodel1, type="II")
## Type II Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## group 0.028479 0.0056958 5 95.373 26.969 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(Rmodel1)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: R.nmol.org.min ~ group + (1 | Run)
## Data: PRdf
##
## REML criterion at convergence: -522
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.8829 -0.4791 -0.0918 0.6491 2.2640
##
## Random effects:
## Groups Name Variance Std.Dev.
## Run (Intercept) 0.0000938 0.009685
## Residual 0.0002112 0.014533
## Number of obs: 103, groups: Run, 4
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -0.002193 0.006108 5.805396 -0.359 0.732326
## groupLarvae3 -0.009453 0.004596 94.121230 -2.057 0.042455 *
## groupLarvae5 -0.017675 0.005047 96.678815 -3.502 0.000701 ***
## groupLarvae6 -0.006828 0.005136 95.446406 -1.329 0.186887
## groupRecruit1 -0.053033 0.005078 96.882613 -10.444 < 2e-16 ***
## groupRecruit2 -0.026840 0.005468 95.310225 -4.908 3.79e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) grpLr3 grpLr5 grpLr6 grpRc1
## groupLarva3 -0.376
## groupLarva5 -0.447 0.455
## groupLarva6 -0.392 0.447 0.478
## groupRecrt1 -0.449 0.453 0.540 0.468
## groupRecrt2 -0.368 0.420 0.449 0.438 0.439
Check assumptions of model for residual normality and variance. Passes assumptions.
qqPlot(residuals(Rmodel1))
## [1] 49 94
leveneTest(residuals(Rmodel1)~group, data=PRdf)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 5 2.0095 0.084 .
## 97
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Conduct post hoc test.
emm<-emmeans(Rmodel1, ~group)
cld(emm)
## group emmean SE df lower.CL upper.CL .group
## Recruit1 -0.05523 0.00593 5.21 -0.0703 -0.04015 1
## Recruit2 -0.02903 0.00656 7.38 -0.0444 -0.01368 2
## Larvae5 -0.01987 0.00593 5.21 -0.0349 -0.00480 23
## Larvae3 -0.01165 0.00617 5.60 -0.0270 0.00371 34
## Larvae6 -0.00902 0.00629 6.24 -0.0243 0.00623 34
## Larvae2 -0.00219 0.00617 5.60 -0.0176 0.01317 4
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 6 estimates
## significance level used: alpha = 0.05
## NOTE: Compact letter displays can be misleading
## because they show NON-findings rather than findings.
## Consider using 'pairs()', 'pwpp()', or 'pwpm()' instead.
Build linear mixed effect model and examine for Photosynthesis
Pmodel1<-lmer(P.nmol.org.min~group + (1|Run), data=PRdf) #run nested within day
anova(Pmodel1, type="II")
## Type II Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## group 0.010495 0.0020989 5 95.411 6.4194 3.496e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(Pmodel1)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: P.nmol.org.min ~ group + (1 | Run)
## Data: PRdf
##
## REML criterion at convergence: -480
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.9914 -0.5376 0.1482 0.5734 2.2312
##
## Random effects:
## Groups Name Variance Std.Dev.
## Run (Intercept) 0.0001243 0.01115
## Residual 0.0003270 0.01808
## Number of obs: 103, groups: Run, 4
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.009563 0.007239 6.048001 1.321 0.234277
## groupLarvae3 0.007303 0.005718 94.053368 1.277 0.204702
## groupLarvae5 0.022404 0.006271 96.825444 3.573 0.000552 ***
## groupLarvae6 0.006786 0.006387 95.578132 1.063 0.290678
## groupRecruit1 -0.006544 0.006307 96.972250 -1.038 0.302046
## groupRecruit2 -0.008268 0.006800 95.425418 -1.216 0.227028
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) grpLr3 grpLr5 grpLr6 grpRc1
## groupLarva3 -0.395
## groupLarva5 -0.467 0.456
## groupLarva6 -0.411 0.448 0.478
## groupRecrt1 -0.470 0.453 0.538 0.468
## groupRecrt2 -0.386 0.420 0.449 0.437 0.439
Check assumptions of model for residual normality and variance. Violates variance assumption, return to this.
qqPlot(residuals(Pmodel1))
## [1] 45 72
leveneTest(residuals(Pmodel1)~group, data=PRdf)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 5 2.888 0.01788 *
## 97
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Conduct post hoc test.
emm<-emmeans(Pmodel1, ~group)
cld(emm)
## group emmean SE df lower.CL upper.CL .group
## Recruit2 0.00130 0.00784 8.12 -0.01673 0.0193 1
## Recruit1 0.00302 0.00702 5.59 -0.01447 0.0205 1
## Larvae2 0.00956 0.00733 5.98 -0.00839 0.0275 1
## Larvae6 0.01635 0.00748 6.76 -0.00147 0.0342 12
## Larvae3 0.01687 0.00733 5.98 -0.00108 0.0348 12
## Larvae5 0.03197 0.00702 5.59 0.01448 0.0495 2
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 6 estimates
## significance level used: alpha = 0.05
## NOTE: Compact letter displays can be misleading
## because they show NON-findings rather than findings.
## Consider using 'pairs()', 'pwpp()', or 'pwpm()' instead.
Build linear mixed effect model and examine for GP
GPmodel1<-lmer(GP.nmol.org.min~group + (1|Run), data=PRdf) #run nested within day
anova(GPmodel1, type="II")
## Type II Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## group 0.022737 0.0045475 5 95.131 7.3014 7.918e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(GPmodel1)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: GP.nmol.org.min ~ group + (1 | Run)
## Data: PRdf
##
## REML criterion at convergence: -416.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.3168 -0.5676 -0.0800 0.6846 3.8487
##
## Random effects:
## Groups Name Variance Std.Dev.
## Run (Intercept) 0.0003902 0.01975
## Residual 0.0006228 0.02496
## Number of obs: 103, groups: Run, 4
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.011687 0.011785 4.970329 0.992 0.3672
## groupLarvae3 0.016756 0.007892 94.047008 2.123 0.0364 *
## groupLarvae5 0.040108 0.008691 96.249656 4.615 1.22e-05 ***
## groupLarvae6 0.013307 0.008830 95.076121 1.507 0.1351
## groupRecruit1 0.046599 0.008749 96.524495 5.326 6.55e-07 ***
## groupRecruit2 0.018265 0.009399 94.965810 1.943 0.0549 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) grpLr3 grpLr5 grpLr6 grpRc1
## groupLarva3 -0.335
## groupLarva5 -0.401 0.454
## groupLarva6 -0.350 0.447 0.479
## groupRecrt1 -0.403 0.451 0.542 0.468
## groupRecrt2 -0.329 0.420 0.450 0.439 0.439
Check assumptions of model for residual normality and variance. Violates variance assumption, return to this.
qqPlot(residuals(GPmodel1))
## [1] 49 94
leveneTest(residuals(GPmodel1)~group, data=PRdf)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 5 2.8511 0.0191 *
## 97
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Conduct post hoc test.
emm<-emmeans(GPmodel1, ~group)
cld(emm)
## group emmean SE df lower.CL upper.CL .group
## Larvae2 0.0117 0.0119 4.91 -0.018969 0.0423 1
## Larvae6 0.0250 0.0120 5.31 -0.005423 0.0554 1
## Larvae3 0.0284 0.0119 4.91 -0.002214 0.0591 12
## Recruit2 0.0300 0.0125 6.09 -0.000441 0.0603 12
## Larvae5 0.0518 0.0115 4.56 0.021347 0.0822 23
## Recruit1 0.0583 0.0115 4.56 0.027839 0.0887 3
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 6 estimates
## significance level used: alpha = 0.05
## NOTE: Compact letter displays can be misleading
## because they show NON-findings rather than findings.
## Consider using 'pairs()', 'pwpp()', or 'pwpm()' instead.
Build linear mixed effect model and examine for GP
PRmodel1<-lmer(ratio~group + (1|Run), data=PRdf) #run nested within day
anova(PRmodel1, type="II")
## Type II Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## group 32.798 6.5596 5 95.277 2.5797 0.03109 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(PRmodel1)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: ratio ~ group + (1 | Run)
## Data: PRdf
##
## REML criterion at convergence: 383.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.3375 -0.4941 -0.1241 0.2552 5.8948
##
## Random effects:
## Groups Name Variance Std.Dev.
## Run (Intercept) 0.02423 0.1557
## Residual 2.54277 1.5946
## Number of obs: 103, groups: Run, 4
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.4873 0.3722 22.4757 6.684 9.12e-07 ***
## groupLarvae3 -0.3534 0.5043 95.4512 -0.701 0.4852
## groupLarvae5 0.2136 0.5232 88.9572 0.408 0.6841
## groupLarvae6 0.1867 0.5478 95.6306 0.341 0.7340
## groupRecruit1 -1.3617 0.5233 88.4715 -2.602 0.0109 *
## groupRecruit2 -0.7642 0.5852 96.2361 -1.306 0.1947
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) grpLr3 grpLr5 grpLr6 grpRc1
## groupLarva3 -0.677
## groupLarva5 -0.680 0.482
## groupLarva6 -0.641 0.460 0.456
## groupRecrt1 -0.680 0.482 0.484 0.455
## groupRecrt2 -0.600 0.431 0.427 0.407 0.426
Check assumptions of model for residual normality and variance. Meets assumptions.
qqPlot(residuals(PRmodel1))
## [1] 14 88
leveneTest(residuals(PRmodel1)~group, data=PRdf)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 5 1.2789 0.2792
## 97
Conduct post hoc test.
emm<-emmeans(PRmodel1, ~group)
cld(emm)
## group emmean SE df lower.CL upper.CL .group
## Recruit1 1.13 0.384 41.4 0.350 1.90 1
## Recruit2 1.72 0.475 50.2 0.769 2.68 12
## Larvae3 2.13 0.403 15.4 1.278 2.99 12
## Larvae2 2.49 0.403 15.4 1.631 3.34 12
## Larvae6 2.67 0.428 38.5 1.808 3.54 12
## Larvae5 2.70 0.384 41.4 1.925 3.48 2
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 6 estimates
## significance level used: alpha = 0.05
## NOTE: Compact letter displays can be misleading
## because they show NON-findings rather than findings.
## Consider using 'pairs()', 'pwpp()', or 'pwpm()' instead.
Generate summary of all respiration data.
summary<-PRdf%>%
group_by(group)%>%
summarise(N=length(R.nmol.org.min),
Mean_R=mean(R.nmol.org.min),
SE_R=sd(R.nmol.org.min)/sqrt(length(R.nmol.org.min)),
Mean_P=mean(P.nmol.org.min),
SE_P=sd(P.nmol.org.min)/sqrt(length(P.nmol.org.min)),
Mean_GP=mean(GP.nmol.org.min),
SE_GP=sd(GP.nmol.org.min)/sqrt(length(GP.nmol.org.min)),
Mean_PR=mean(ratio),
SE_PR=sd(ratio)/sqrt(length(ratio)))
summary%>%
write_csv(., "Mcap2020/Output/Respiration/mean_respiration.csv")
Generate plots by lifestage and days post fertilization.
Generate dot plot.
Rplot1<-ggplot(data=PRdata, aes(x=dpf, y=R.nmol.org.min, colour=Lifestage, group=interaction(Lifestage, dpf))) +
geom_boxplot(position=position_dodge(0.9), lwd=1)+
geom_jitter(shape=16, position=position_dodge(0.9), size=3)+
geom_hline(yintercept=0, linetype="dashed", color="black", size=0.75)+
scale_color_brewer(palette = "Dark2", breaks=c("Larvae", "Recruit", "Juvenile"))+
scale_y_continuous(limits=c(-0.1, 0.1), labels = scales::number_format(accuracy = 0.01, decimal.mark = '.'))+
xlab("Days post fertilization")+
ylab(expression(bold(paste("R (nmol ", O[2], " larva"^-1, "min"^-1, ")")))) +
theme_classic() +
theme(text = element_text(size = 18, color="black"),
axis.text = element_text(size = 18, color="black"),
legend.position = "top",
axis.title = element_text(size = 18, color="black", face="bold"),
legend.title=element_blank(),
legend.text = element_text(size=18),
plot.margin=unit(c(1,1,1,1), "cm"),
axis.title.y = element_text(margin = margin(t = 0, r = 1, b = 0, l = 1)),
axis.title.x = element_text(margin = margin(t = 3, r = 0, b = 0, l = 0)));Rplot1
Rplot2<-Rplot1+theme(legend.position="none")
#ggsave("Mcap2020/Figures/Respiration/mcap_R.pdf", plot=Rplot1, height=8, width=7, units = c("in"), dpi=300) #output figure
Generate mean table
meanR <- plyr::ddply(PRdata, c("Lifestage", "dpf"), summarise,
N = length(R.nmol.org.min[!is.na(R.nmol.org.min)]),
mean = mean(R.nmol.org.min, na.rm=TRUE),
sd = sd(R.nmol.org.min, na.rm=TRUE),
se = sd / sqrt(N)
); meanR
## Lifestage dpf N mean sd se
## 1 Embryo 3 20 -0.004191302 0.017736649 0.0039660353
## 2 Juvenile 40 20 -0.045210023 0.019598294 0.0043823118
## 3 Larvae 4 20 -0.013644439 0.004434086 0.0009914917
## 4 Larvae 8 18 -0.019675601 0.005887039 0.0013875884
## 5 Larvae 10 15 -0.011576113 0.018175545 0.0046929056
## 6 Recruit 8 18 -0.055417500 0.027867296 0.0065683846
## 7 Recruit 10 12 -0.031588057 0.017438224 0.0050339817
Pplot1<-ggplot(data=PRdata, aes(x=dpf, y=P.nmol.org.min, colour=Lifestage, group=interaction(Lifestage, dpf))) +
geom_boxplot(position=position_dodge(0.9), lwd=1)+
geom_jitter(shape=16, position=position_dodge(0.9), size=3)+
geom_hline(yintercept=0, linetype="dashed", color="black", size=0.75)+
scale_color_brewer(palette = "Dark2", breaks=c("Larvae", "Recruit", "Juvenile"))+
scale_y_continuous(limits=c(-0.1, 0.3), labels = scales::number_format(accuracy = 0.01, decimal.mark = '.'))+
xlab("Days post fertilization")+
ylab(expression(bold(paste("Net P (nmol ", O[2], " larva"^-1, "min"^-1, ")")))) +
theme_classic() +
theme(text = element_text(size = 18, color="black"),
axis.text = element_text(size = 18, color="black"),
legend.position = "right",
axis.title = element_text(size = 18, color="black", face="bold"),
legend.title=element_blank(),
legend.text = element_text(size=18),
plot.margin=unit(c(1,1,1,1), "cm"),
axis.title.y = element_text(margin = margin(t = 0, r = 1, b = 0, l = 1)),
axis.title.x = element_text(margin = margin(t = 3, r = 0, b = 0, l = 0))); Pplot1
Pplot2<-Pplot1+theme(legend.position="none")
#ggsave("Mcap2020/Figures/Respiration/mcap_P.pdf", plot=Pplot2, height=8, width=7, units = c("in"), dpi=300) #output figure
Generate mean table.
meanP <- plyr::ddply(PRdata, c("Lifestage", "dpf"), summarise,
N = length(P.nmol.org.min[!is.na(P.nmol.org.min)]),
mean = mean(P.nmol.org.min, na.rm=TRUE),
sd = sd(P.nmol.org.min, na.rm=TRUE),
se = sd / sqrt(N)
); meanP
## Lifestage dpf N mean sd se
## 1 Embryo 3 20 0.006335109 0.01133371 0.002534295
## 2 Juvenile 40 20 0.079503148 0.03583541 0.008013042
## 3 Larvae 4 20 0.013637833 0.01476790 0.003302203
## 4 Larvae 8 18 0.032156372 0.01176212 0.002772357
## 5 Larvae 10 15 0.018565242 0.02653452 0.006851184
## 6 Recruit 8 18 0.002830864 0.02878711 0.006785188
## 7 Recruit 10 12 0.003511298 0.02449710 0.007071705
GPplot1<-ggplot(data=PRdata, aes(x=dpf, y=GP.nmol.org.min, colour=Lifestage, group=interaction(Lifestage, dpf))) +
geom_boxplot(position=position_dodge(0.9), lwd=1)+
geom_jitter(shape=16, position=position_dodge(0.9), size=3)+
geom_hline(yintercept=0, linetype="dashed", color="black", size=0.75)+
scale_color_brewer(palette = "Dark2", breaks=c("Larvae", "Recruit", "Juvenile"))+
scale_y_continuous(limits=c(-0.1, 0.3), labels = scales::number_format(accuracy = 0.01, decimal.mark = '.'))+
xlab("Days post fertilization")+
ylab(expression(bold(paste("Gross P (nmol ", O[2], " larva"^-1, "min"^-1, ")")))) +
theme_classic() +
theme(text = element_text(size = 18, color="black"),
axis.text = element_text(size = 18, color="black"),
legend.position = "right",
axis.title = element_text(size = 18, color="black", face="bold"),
legend.title=element_blank(),
legend.text = element_text(size=18),
plot.margin=unit(c(1,1,1,1), "cm"),
axis.title.y = element_text(margin = margin(t = 0, r = 1, b = 0, l = 1)),
axis.title.x = element_text(margin = margin(t = 3, r = 0, b = 0, l = 0))); GPplot1
GPplot2<-GPplot1+theme(legend.position="none")
#ggsave("Mcap2020/Figures/Respiration/mcap_GP.pdf", plot=GPplot2, height=8, width=7, units = c("in"), dpi=300) #output figure
Generate mean table.
meanGP <- plyr::ddply(PRdata, c("Lifestage", "dpf"), summarise,
N = length(GP.nmol.org.min[!is.na(GP.nmol.org.min)]),
mean = mean(GP.nmol.org.min, na.rm=TRUE),
sd = sd(GP.nmol.org.min, na.rm=TRUE),
se = sd / sqrt(N)
); meanGP
## Lifestage dpf N mean sd se
## 1 Embryo 3 20 0.01052641 0.02591563 0.005794912
## 2 Juvenile 40 20 0.12471317 0.04749814 0.010620907
## 3 Larvae 4 20 0.02728227 0.01561791 0.003492270
## 4 Larvae 8 18 0.05183197 0.01398628 0.003296597
## 5 Larvae 10 15 0.03014136 0.04341109 0.011208695
## 6 Recruit 8 18 0.05824836 0.03865966 0.009112168
## 7 Recruit 10 12 0.03509936 0.03818919 0.011024270
PRplot1<-ggplot(data=PRdata, aes(x=dpf, y=ratio, colour=Lifestage, group=interaction(Lifestage, dpf))) +
geom_boxplot(position=position_dodge(0.9), lwd=1)+
geom_jitter(shape=16, position=position_dodge(0.9), size=3)+
scale_color_brewer(palette = "Dark2", breaks=c("Larvae", "Recruit", "Juvenile"))+
geom_hline(yintercept=1, linetype="dashed", color="black", size=0.75)+
scale_y_continuous(limits=c(-2, 6), labels = scales::number_format(accuracy = 0.01, decimal.mark = '.'))+
ylab(expression(bold(paste("P : R ")))) +
xlab("Days Post Fertilization") +
theme_classic() +
theme(text = element_text(size = 18, color="black"),
axis.text = element_text(size = 18, color="black"),
legend.position = "top",
axis.title = element_text(size = 18, color="black", face="bold"),
legend.title=element_blank(),
legend.text = element_text(size=22),
plot.margin=unit(c(1,1,1,1), "cm"),
axis.title.y = element_text(margin = margin(t = 0, r = 1, b = 0, l = 1)),
axis.title.x = element_text(margin = margin(t = 3, r = 0, b = 0, l = 0))); PRplot1
PRplot2<-PRplot1+theme(legend.position="right")
#ggsave("Mcap2020/Figures/Respiration/mcap_PR.pdf", plot=PRplot1, height=8, width=8, units = c("in"), dpi=300) #output figure
Generate mean table.
meanPR <- plyr::ddply(PRdata, c("Lifestage", "dpf"), summarise,
N = length(ratio[!is.na(ratio)]),
mean = mean(ratio, na.rm=TRUE),
sd = sd(ratio, na.rm=TRUE),
se = sd / sqrt(N)
); meanPR
## Lifestage dpf N mean sd se
## 1 Embryo 3 20 2.492866 2.6746501 0.5980699
## 2 Juvenile 40 20 3.002808 1.0719306 0.2396910
## 3 Larvae 4 20 2.139497 1.0377752 0.2320536
## 4 Larvae 8 18 2.705802 0.7656908 0.1804751
## 5 Larvae 10 15 2.680163 1.9589490 0.5057985
## 6 Recruit 8 18 1.120693 0.6434859 0.1516711
## 7 Recruit 10 12 1.729208 1.3888877 0.4009374
Combine plots above into a single panel.
mcap_respiration<-plot_grid(Rplot2, Pplot2, GPplot2, PRplot2, labels = c("A", "B", "C", "D"), label_size=18, ncol=4, nrow=1, rel_heights= c(1,1,1,1), rel_widths = c(0.75,0.75,0.75,1), align="h")
#ggsave(filename="Mcap2020/Figures/Respiration/mcap_PR_panel.pdf", plot=mcap_respiration, dpi=300, width=24, height=6, units="in")
Generate a plot for fused corals separately to examine respirometry normalized to the number polyps fused together within a colony.
FuseR1<-ggplot(data=fused_oxygen, aes(x=as.factor(Org.Number), y=R.nmol.org.min, group=interaction(Org.Number))) +
geom_boxplot(position=position_dodge(0.9), lwd=1)+
geom_jitter(shape=16, position=position_dodge(0.9), size=3)+
geom_hline(yintercept=0, linetype="dashed", color="black", size=0.75)+
scale_y_continuous(limits=c(-0.1, 0.01), labels = scales::number_format(accuracy = 0.01, decimal.mark = '.'))+
xlab("Number of Fused Polyps")+
ylab(expression(bold(paste("R (nmol ", O[2], " polyp"^-1, "min"^-1, ")")))) +
theme_classic() +
theme(text = element_text(size = 18, color="black"),
axis.text = element_text(size = 18, color="black"),
legend.position = "right",
axis.title = element_text(size = 18, color="black", face="bold"),
legend.title=element_blank(),
legend.text = element_text(size=18),
plot.margin=unit(c(1,1,1,1), "cm"),
axis.title.y = element_text(margin = margin(t = 0, r = 1, b = 0, l = 1)),
axis.title.x = element_text(margin = margin(t = 3, r = 0, b = 0, l = 0)));FuseR1
FuseP1<-ggplot(data=fused_oxygen, aes(x=as.factor(Org.Number), y=P.nmol.org.min, group=interaction(Org.Number))) +
geom_boxplot(position=position_dodge(0.9), lwd=1)+
geom_jitter(shape=16, position=position_dodge(0.9), size=3)+
geom_hline(yintercept=0, linetype="dashed", color="black", size=0.75)+
scale_y_continuous(limits=c(-0.1, 0.4), labels = scales::number_format(accuracy = 0.01, decimal.mark = '.'))+
xlab("Number of Fused Polyps")+
ylab(expression(bold(paste("Net P (nmol ", O[2], " polyp"^-1, "min"^-1, ")")))) +
theme_classic() +
theme(text = element_text(size = 18, color="black"),
axis.text = element_text(size = 18, color="black"),
legend.position = "right",
axis.title = element_text(size = 18, color="black", face="bold"),
legend.title=element_blank(),
legend.text = element_text(size=18),
plot.margin=unit(c(1,1,1,1), "cm"),
axis.title.y = element_text(margin = margin(t = 0, r = 1, b = 0, l = 1)),
axis.title.x = element_text(margin = margin(t = 3, r = 0, b = 0, l = 0)));FuseP1
FuseGP1<-ggplot(data=fused_oxygen, aes(x=as.factor(Org.Number), y=GP.nmol.org.min, group=interaction(Org.Number))) +
geom_boxplot(position=position_dodge(0.9), lwd=1)+
geom_jitter(shape=16, position=position_dodge(0.9), size=3)+
geom_hline(yintercept=0, linetype="dashed", color="black", size=0.75)+
scale_y_continuous(limits=c(-0.1, 0.4), labels = scales::number_format(accuracy = 0.01, decimal.mark = '.'))+
xlab("Number of Fused Polyps")+
ylab(expression(bold(paste("Gross P (nmol ", O[2], " polyp"^-1, "min"^-1, ")")))) +
theme_classic() +
theme(text = element_text(size = 18, color="black"),
axis.text = element_text(size = 18, color="black"),
legend.position = "right",
axis.title = element_text(size = 18, color="black", face="bold"),
legend.title=element_blank(),
legend.text = element_text(size=18),
plot.margin=unit(c(1,1,1,1), "cm"),
axis.title.y = element_text(margin = margin(t = 0, r = 1, b = 0, l = 1)),
axis.title.x = element_text(margin = margin(t = 3, r = 0, b = 0, l = 0)));FuseGP1
FusePR1<-ggplot(data=fused_oxygen, aes(x=as.factor(Org.Number), y=ratio, group=interaction(Org.Number))) +
geom_boxplot(position=position_dodge(0.9), lwd=1)+
geom_jitter(shape=16, position=position_dodge(0.9), size=3)+
geom_hline(yintercept=1, linetype="dashed", color="black", size=0.75)+
scale_y_continuous(limits=c(-1, 8), labels = scales::number_format(accuracy = 0.01, decimal.mark = '.'))+
xlab("Number of Fused Polyps")+
ylab(expression(bold(paste("P : R")))) +
theme_classic() +
theme(text = element_text(size = 18, color="black"),
axis.text = element_text(size = 18, color="black"),
legend.position = "right",
axis.title = element_text(size = 18, color="black", face="bold"),
legend.title=element_blank(),
legend.text = element_text(size=18),
plot.margin=unit(c(1,1,1,1), "cm"),
axis.title.y = element_text(margin = margin(t = 0, r = 1, b = 0, l = 1)),
axis.title.x = element_text(margin = margin(t = 3, r = 0, b = 0, l = 0)));FusePR1
Fuse_c_R1<-ggplot(data=fused_oxygen, aes(x=as.factor(Org.Number), y=(R.nmol.org.min*Org.Number), group=interaction(Org.Number))) +
geom_boxplot(position=position_dodge(0.9), lwd=1)+
geom_jitter(shape=16, position=position_dodge(0.9), size=3)+
geom_hline(yintercept=0, linetype="dashed", color="black", size=0.75)+
scale_y_continuous(limits=c(-0.4, 0.1), labels = scales::number_format(accuracy = 0.01, decimal.mark = '.'))+
xlab("Number of Fused Polyps")+
ylab(expression(bold(paste("R (nmol ", O[2], " polyp"^-1, "min"^-1, ")")))) +
theme_classic() +
theme(text = element_text(size = 18, color="black"),
axis.text = element_text(size = 18, color="black"),
legend.position = "right",
axis.title = element_text(size = 18, color="black", face="bold"),
legend.title=element_blank(),
legend.text = element_text(size=18),
plot.margin=unit(c(1,1,1,1), "cm"),
axis.title.y = element_text(margin = margin(t = 0, r = 1, b = 0, l = 1)),
axis.title.x = element_text(margin = margin(t = 3, r = 0, b = 0, l = 0)));Fuse_c_R1
Fuse_c_P1<-ggplot(data=fused_oxygen, aes(x=as.factor(Org.Number), y=(P.nmol.org.min*Org.Number), group=interaction(Org.Number))) +
geom_boxplot(position=position_dodge(0.9), lwd=1)+
geom_jitter(shape=16, position=position_dodge(0.9), size=3)+
geom_hline(yintercept=0, linetype="dashed", color="black", size=0.75)+
scale_y_continuous(limits=c(-0.1, 0.8), labels = scales::number_format(accuracy = 0.01, decimal.mark = '.'))+
xlab("Number of Fused Polyps")+
ylab(expression(bold(paste("Net P (nmol ", O[2], " polyp"^-1, "min"^-1, ")")))) +
theme_classic() +
theme(text = element_text(size = 18, color="black"),
axis.text = element_text(size = 18, color="black"),
legend.position = "right",
axis.title = element_text(size = 18, color="black", face="bold"),
legend.title=element_blank(),
legend.text = element_text(size=18),
plot.margin=unit(c(1,1,1,1), "cm"),
axis.title.y = element_text(margin = margin(t = 0, r = 1, b = 0, l = 1)),
axis.title.x = element_text(margin = margin(t = 3, r = 0, b = 0, l = 0)));Fuse_c_P1
Fuse_c_GP1<-ggplot(data=fused_oxygen, aes(x=as.factor(Org.Number), y=(GP.nmol.org.min*Org.Number), group=interaction(Org.Number))) +
geom_boxplot(position=position_dodge(0.9), lwd=1)+
geom_jitter(shape=16, position=position_dodge(0.9), size=3)+
geom_hline(yintercept=0, linetype="dashed", color="black", size=0.75)+
scale_y_continuous(limits=c(-0.1, 0.8), labels = scales::number_format(accuracy = 0.01, decimal.mark = '.'))+
xlab("Number of Fused Polyps")+
ylab(expression(bold(paste("Gross P (nmol ", O[2], " polyp"^-1, "min"^-1, ")")))) +
theme_classic() +
theme(text = element_text(size = 18, color="black"),
axis.text = element_text(size = 18, color="black"),
legend.position = "right",
axis.title = element_text(size = 18, color="black", face="bold"),
legend.title=element_blank(),
legend.text = element_text(size=18),
plot.margin=unit(c(1,1,1,1), "cm"),
axis.title.y = element_text(margin = margin(t = 0, r = 1, b = 0, l = 1)),
axis.title.x = element_text(margin = margin(t = 3, r = 0, b = 0, l = 0)));Fuse_c_GP1
Fuse_c_PR1<-ggplot(data=fused_oxygen, aes(x=as.factor(Org.Number), y=(abs(fused_oxygen$GP.nmol.org.min*Org.Number)/abs(fused_oxygen$R.nmol.org.min*Org.Number)), group=interaction(Org.Number))) +
geom_boxplot(position=position_dodge(0.9), lwd=1)+
geom_jitter(shape=16, position=position_dodge(0.9), size=3)+
geom_hline(yintercept=1, linetype="dashed", color="black", size=0.75)+
scale_y_continuous(limits=c(-1, 10), labels = scales::number_format(accuracy = 0.01, decimal.mark = '.'))+
xlab("Number of Fused Polyps")+
ylab(expression(bold(paste("P : R")))) +
theme_classic() +
theme(text = element_text(size = 18, color="black"),
axis.text = element_text(size = 18, color="black"),
legend.position = "right",
axis.title = element_text(size = 18, color="black", face="bold"),
legend.title=element_blank(),
legend.text = element_text(size=18),
plot.margin=unit(c(1,1,1,1), "cm"),
axis.title.y = element_text(margin = margin(t = 0, r = 1, b = 0, l = 1)),
axis.title.x = element_text(margin = margin(t = 3, r = 0, b = 0, l = 0)));Fuse_c_PR1
Combine plots above into a single panel.
mcap_fused_respiration<-plot_grid(FuseR1, FuseP1, FuseGP1, FusePR1, Fuse_c_R1, Fuse_c_P1, Fuse_c_GP1, Fuse_c_PR1, labels = c("Polyp - Level Responses", "", "", "", "Colony - Level Responses", "", "", ""), label_size = 18, ncol=4, nrow=2, rel_heights= c(1), rel_widths = c(1), align="h")
ggsave(filename="Mcap2020/Figures/Respiration/mcap_fused_panel.pdf", plot=mcap_fused_respiration, dpi=300, width=24, height=12, units="in")
Make a summary table to displays means and standard errors for plotting.
colony_v_polypsR <- plyr::ddply(fused_oxygen, c("Org.Number"), summarise,
N_r_polyps = length(R.nmol.org.min[!is.na(R.nmol.org.min)]),
mean_r_polyps = mean(R.nmol.org.min, na.rm=TRUE),
sd_r_polyps = sd(R.nmol.org.min, na.rm=TRUE),
se_r_polyps = sd_r_polyps / sqrt(N_r_polyps),
N_r_colony = length(R.nmol.org.min[!is.na(R.nmol.org.min)]),
mean_r_colony = mean((Org.Number*R.nmol.org.min), na.rm=TRUE),
sd_r_colony = sd(Org.Number*R.nmol.org.min, na.rm=TRUE),
se_r_colony = sd_r_colony / sqrt(N_r_colony)
)
Plot a line plot with two axis. Left axis: polyp specific respiration. Right axis: colony total respiration.
R_c_vs_p<-ggplot(data=colony_v_polypsR, aes(x=as.factor(Org.Number), group=interaction(Org.Number))) +
geom_point(aes(y=mean_r_polyps), size=3, position=position_dodge(0.01), color="orange") + #polyp
geom_errorbar(aes(ymin=mean_r_polyps-se_r_polyps, ymax=mean_r_polyps+se_r_polyps), width=0, linetype="solid", color="orange", position=position_dodge(0.01), size=0.5) + #polyp
geom_point(aes(y=mean_r_colony), size=3, position=position_dodge(0.06), color="purple") + #polyp
geom_errorbar(aes(ymin=mean_r_colony-se_r_colony, ymax=mean_r_colony+se_r_colony), width=0, linetype="solid", color="purple", position=position_dodge(0.06), size=0.5) + #polyp
geom_hline(yintercept=0, linetype="dashed", color="black", size=0.75)+
scale_y_continuous(
name = "Polyp Respiration", # Features of the first axis
sec.axis = sec_axis(~.*1, name="Colony Total Respiration"), # Add a second axis and specify its features
limits=c(-0.4, 0.1),
labels = scales::number_format(accuracy = 0.01, decimal.mark = '.')
) +
xlab("Number of Fused Polyps")+
#ylab(expression(bold(paste("Respiration: nmol polyp"^-1, "min"^-1)))) +
theme_classic() +
theme(text = element_text(size = 18, color="black"),
axis.text = element_text(size = 18, color="black"),
legend.position = "right",
axis.title = element_text(size = 18, color="black", face="bold"),
legend.title=element_blank(),
legend.text = element_text(size=18),
plot.margin=unit(c(1,1,1,1), "cm"),
axis.title.y = element_text(margin = margin(t = 0, r = 1, b = 0, l = 1), color="orange", face="bold"),
axis.title.y.right = element_text(margin = margin(t = 0, r = 1, b = 1, l = 2), color="purple", face="bold"),
axis.title.x = element_text(margin = margin(t = 3, r = 0, b = 0, l = 0)));R_c_vs_p
Make a summary table to displays means and standard errors for plotting.
colony_v_polypsP <- plyr::ddply(fused_oxygen, c("Org.Number"), summarise,
N_p_polyps = length(P.nmol.org.min[!is.na(P.nmol.org.min)]),
mean_p_polyps = mean(P.nmol.org.min, na.rm=TRUE),
sd_p_polyps = sd(P.nmol.org.min, na.rm=TRUE),
se_p_polyps = sd_p_polyps / sqrt(N_p_polyps),
N_p_colony = length(P.nmol.org.min[!is.na(P.nmol.org.min)]),
mean_p_colony = mean((Org.Number*P.nmol.org.min), na.rm=TRUE),
sd_p_colony = sd(Org.Number*P.nmol.org.min, na.rm=TRUE),
se_p_colony = sd_p_colony / sqrt(N_p_colony)
)
Plot a line plot with two axis. Left axis: polyp specific net photosynthesis Right axis: colony total net photosynthesis.
P_c_vs_p<-ggplot(data=colony_v_polypsP, aes(x=as.factor(Org.Number), group=interaction(Org.Number))) +
geom_point(aes(y=mean_p_polyps), size=3, position=position_dodge(0.01), color="orange") + #polyp
geom_errorbar(aes(ymin=mean_p_polyps-se_p_polyps, ymax=mean_p_polyps+se_p_polyps), width=0, linetype="solid", color="orange", position=position_dodge(0.01), size=0.5) + #polyp
geom_point(aes(y=mean_p_colony), size=3, position=position_dodge(0.06), color="purple") + #polyp
geom_errorbar(aes(ymin=mean_p_colony-se_p_colony, ymax=mean_p_colony+se_p_colony), width=0, linetype="solid", color="purple", position=position_dodge(0.06), size=0.5) + #polyp
geom_hline(yintercept=0, linetype="dashed", color="black", size=0.75)+
scale_y_continuous(
name = "Polyp Net Photosynthesis", # Features of the first axis
sec.axis = sec_axis(~.*1, name="Colony Total Net Photosynthesis"), # Add a second axis and specify its features
limits=c(-0.1, 0.8),
labels = scales::number_format(accuracy = 0.01, decimal.mark = '.')
) +
xlab("Number of Fused Polyps")+
#ylab(expression(bold(paste("Respiration: nmol polyp"^-1, "min"^-1)))) +
theme_classic() +
theme(text = element_text(size = 18, color="black"),
axis.text = element_text(size = 18, color="black"),
legend.position = "right",
axis.title = element_text(size = 18, color="black", face="bold"),
legend.title=element_blank(),
legend.text = element_text(size=18),
plot.margin=unit(c(1,1,1,1), "cm"),
axis.title.y = element_text(margin = margin(t = 0, r = 1, b = 0, l = 1), color="orange", face="bold"),
axis.title.y.right = element_text(margin = margin(t = 0, r = 1, b = 1, l = 2), color="purple", face="bold"),
axis.title.x = element_text(margin = margin(t = 3, r = 0, b = 0, l = 0)));P_c_vs_p
Make a summary table to displays means and standard errors for plotting.
colony_v_polypsGP <- plyr::ddply(fused_oxygen, c("Org.Number"), summarise,
N_gp_polyps = length(GP.nmol.org.min[!is.na(GP.nmol.org.min)]),
mean_gp_polyps = mean(GP.nmol.org.min, na.rm=TRUE),
sd_gp_polyps = sd(GP.nmol.org.min, na.rm=TRUE),
se_gp_polyps = sd_gp_polyps / sqrt(N_gp_polyps),
N_gp_colony = length(GP.nmol.org.min[!is.na(GP.nmol.org.min)]),
mean_gp_colony = mean((Org.Number*GP.nmol.org.min), na.rm=TRUE),
sd_gp_colony = sd(Org.Number*GP.nmol.org.min, na.rm=TRUE),
se_gp_colony = sd_gp_colony / sqrt(N_gp_colony)
)
Plot a line plot with two axis. Left axis: polyp specific net photosynthesis Right axis: colony total net photosynthesis.
GP_c_vs_p<-ggplot(data=colony_v_polypsGP, aes(x=as.factor(Org.Number), group=interaction(Org.Number))) +
geom_point(aes(y=mean_gp_polyps), size=3, position=position_dodge(0.01), color="orange") + #polyp
geom_errorbar(aes(ymin=mean_gp_polyps-se_gp_polyps, ymax=mean_gp_polyps+se_gp_polyps), width=0, linetype="solid", color="orange", position=position_dodge(0.01), size=0.5) + #polyp
geom_point(aes(y=mean_gp_colony), size=3, position=position_dodge(0.06), color="purple") + #polyp
geom_errorbar(aes(ymin=mean_gp_colony-se_gp_colony, ymax=mean_gp_colony+se_gp_colony), width=0, linetype="solid", color="purple", position=position_dodge(0.06), size=0.5) + #polyp
geom_hline(yintercept=0, linetype="dashed", color="black", size=0.75)+
scale_y_continuous(
name = "Polyp Gross Photosynthesis", # Features of the first axis
sec.axis = sec_axis(~.*1, name="Colony Total Gross Photosynthesis"), # Add a second axis and specify its features
limits=c(-0.1, 0.8),
labels = scales::number_format(accuracy = 0.01, decimal.mark = '.')
) +
xlab("Number of Fused Polyps")+
#ylab(expression(bold(paste("Respiration: nmol polyp"^-1, "min"^-1)))) +
theme_classic() +
theme(text = element_text(size = 18, color="black"),
axis.text = element_text(size = 18, color="black"),
legend.position = "right",
axis.title = element_text(size = 18, color="black", face="bold"),
legend.title=element_blank(),
legend.text = element_text(size=18),
plot.margin=unit(c(1,1,1,1), "cm"),
axis.title.y = element_text(margin = margin(t = 0, r = 1, b = 0, l = 1), color="orange", face="bold"),
axis.title.y.right = element_text(margin = margin(t = 0, r = 1, b = 1, l = 2), color="purple", face="bold"),
axis.title.x = element_text(margin = margin(t = 3, r = 0, b = 0, l = 0)));GP_c_vs_p
P:R is the same when calculated for polyp or colony rates. Display P:R means in plot.
Generate mean plot.
meanPR_fuse <- plyr::ddply(fused_oxygen, c("Org.Number"), summarise,
N = length(ratio[!is.na(ratio)]),
mean = mean(ratio, na.rm=TRUE),
sd = sd(ratio, na.rm=TRUE),
se = sd / sqrt(N)
)
PRplot_fuse<-ggplot(data=meanPR_fuse, aes(x=as.factor(Org.Number), y=mean, group=interaction(Org.Number))) +
geom_point(color="darkgray", size=3, position=position_dodge(0.3)) +
geom_errorbar(aes(ymin=mean-se, ymax=mean+se), width=0, linetype="solid", position=position_dodge(0.3), size=1) +
geom_hline(yintercept=1, linetype="dashed", color="black", size=0.75)+
scale_y_continuous(limits=c(-1, 8), labels = scales::number_format(accuracy = 1, decimal.mark = '.'))+
xlab("Number of Fused Polyps")+
ylab(expression(bold(paste("Gross Photosynthesis: Respiration")))) +
theme_classic() +
theme(text = element_text(size = 18, color="black"),
axis.text = element_text(size = 18, color="black"),
legend.position = "right",
axis.title = element_text(size = 18, color="black", face="bold"),
legend.title=element_blank(),
legend.text = element_text(size=18),
plot.margin=unit(c(1,1,1,1), "cm"),
axis.title.y = element_text(margin = margin(t = 0, r = 1, b = 0, l = 1)),
axis.title.x = element_text(margin = margin(t = 3, r = 0, b = 0, l = 0)));PRplot_fuse
Combine plots above into a single panel.
mcap_fused_colony_vs_polyp<-plot_grid(R_c_vs_p, P_c_vs_p, GP_c_vs_p, PRplot_fuse, labels = c("Respiration", "Net Photosynthesis", "Gross Photosynthesis", "P:R"), label_size = 16, ncol=4, nrow=1, rel_heights= c(1,1,1,1), rel_widths = c(1,1,1,0.8), align="h")
ggsave(filename="Mcap2020/Figures/Respiration/mcap_fused_colony_vs_polyp_panel.pdf", plot=mcap_fused_colony_vs_polyp, dpi=300, width=24, height=6, units="in")
Build linear mixed effect model and examine for Respiration.
Rmodel1<-lmer(R.nmol.org.min~Lifestage*dpf + (1|Run/dpf), data=PRdata) #run nested within day
anova(Rmodel1, type="II")
## Type II Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Lifestage 0.0149887 0.0049962 3 8.786 27.1566 8.829e-05 ***
## dpf 0.0007022 0.0003511 2 6.163 1.9085 0.2264
## Lifestage:dpf 0.0008865 0.0008865 1 108.133 4.8184 0.0303 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Check assumptions of model for residual normality and variance. Variance is violated, revisit this assumption.
qqPlot(residuals(Rmodel1))
## [1] 69 38
hist(residuals(Rmodel1))
leveneTest(residuals(Rmodel1)~Lifestage * dpf, data=PRdata)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 6 3.6135 0.002538 **
## 116
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
For respiration, lifestage and lifestage:days post fertilization are significant. Respiration increased over larval development and was higher for recruits at day 8. This is the closest timepoint to metamorphosis, which may suggest that respiratory demand is high immediately after settlement and then reduces by day 10. Juveniles at 40 days post settlement had high respiration, indicating higher demand with larger size.
Build linear mixed effect model and examine for Net Photosynthesis.
Pmodel1<-lmer(P.nmol.org.min~Lifestage*dpf + (1|Run/dpf), data=PRdata)
anova(Pmodel1, type="II")
## Type II Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Lifestage 0.034972 0.0116572 3 7.806 25.5633 0.0002152 ***
## dpf 0.001209 0.0006044 2 5.949 1.3254 0.3341582
## Lifestage:dpf 0.000709 0.0007094 1 108.287 1.5556 0.2149945
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Check assumptions of model for residual normality and variance. Variance is just barely violated, revisit this assumption.
qqPlot(residuals(Pmodel1))
## [1] 21 27
hist(residuals(Pmodel1))
leveneTest(residuals(Pmodel1)~Lifestage * dpf, data=PRdata)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 6 4.7333 0.0002394 ***
## 116
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
For net photosynthesis, lifestage is significant. Photosynthesis appears to be higher in larvae (esp. at day 8), but reduces in early recruits (day 8 and day 10), which may suggest that symbiosis is “upset” or “reset” after metamorphosis and takes time to build to high levels seen in 40 day old juveniles. This supports our hypotheses that symbiosis integration occurs after settlement during early juvenile development.
Build linear mixed effect model and examine for Gross Photosynthesis.
GPmodel1<-lmer(GP.nmol.org.min~Lifestage*dpf + (1|Run/dpf), data=PRdata)
anova(GPmodel1, type="II")
## Type II Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Lifestage 0.0244856 0.0081619 3 8.158 11.1365 0.00297 **
## dpf 0.0027942 0.0013971 2 5.942 1.9063 0.22931
## Lifestage:dpf 0.0000147 0.0000147 1 108.160 0.0201 0.88762
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Check assumptions of model for residual normality and variance.
qqPlot(residuals(GPmodel1))
## [1] 38 69
hist(residuals(GPmodel1))
leveneTest(residuals(GPmodel1)~Lifestage * dpf, data=PRdata)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 6 4.1392 0.0008364 ***
## 116
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
For gross photosynthesis, lifestage is significant. Gross photosynthesis does not appear to reduce as much as net photosynthesis after settlement, so this may indicate that the heavy respiratory demand is what is leading to decreased net photosynthetic output and therefore less available for growth/calcification at this life stage. Gross photosynthesis is still highest in juveniles, which supports that symbiosis is active and integrated in older juveniles.
Build linear mixed effect model and examine for P:R ratio (gross photosynthesis : respiration)
PRmodel1<-lmer(ratio~Lifestage*dpf + (1|Run/dpf), data=PRdata)
anova(PRmodel1, type="II")
## Type II Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Lifestage 33.230 11.0767 3 57.576 4.7804 0.004827 **
## dpf 4.938 2.4691 2 81.324 1.0656 0.349280
## Lifestage:dpf 1.508 1.5075 1 114.313 0.6506 0.421567
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Check assumptions of model for residual normality and variance. Assumptions met.
qqPlot(residuals(PRmodel1))
## [1] 14 108
hist(residuals(PRmodel1))
leveneTest(residuals(PRmodel1)~Lifestage * dpf, data=PRdata)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 6 1.1971 0.3128
## 116
Lifestage is significant in analysis of P:R. P:R values are barely above 1 for larvae and are lowest for recruits, emphasizing the high respiratory demand not met by photosynthesis at this stage. Ratio is higher in juveniles, providing more energy for growth and calcification.
Analyze respiration at the level of the whole colony, or normalized at polyp-specific rates within each colony. Analyze across total number of polyps fused together.
Build linear mixed effect model and examine for Respiration at the colony level.
Rmodel2<-lmer((R.nmol.org.min*Org.Number)~Org.Number + (1|Run/dpf), data=fused_oxygen) #run nested within day
anova(Rmodel2, type="II")
## Type II Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Org.Number 0.037674 0.037674 1 18 6.4345 0.02067 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Check assumptions of model for residual normality and variance. Assumptions are met.
qqPlot(residuals(Rmodel2))
## 38 30
## 18 10
hist(residuals(Rmodel2))
leveneTest(residuals(Rmodel2)~as.factor(Org.Number), data=fused_oxygen)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 7 2.3049 0.0975 .
## 12
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Build linear mixed effect model and examine for Respiration at the polyp level.
Rmodel3<-lmer(R.nmol.org.min~Org.Number + (1|Run/dpf), data=fused_oxygen) #run nested within day
anova(Rmodel3, type="II")
## Type II Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Org.Number 0.0017085 0.0017085 1 18 5.5023 0.03065 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Check assumptions of model for residual normality and variance. Assumptions met.
qqPlot(residuals(Rmodel3))
## 38 30
## 18 10
hist(residuals(Rmodel3))
leveneTest(residuals(Rmodel3)~as.factor(Org.Number), data=fused_oxygen)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 7 1.0323 0.4579
## 12
Display linear regression of 1) colony respiration and 2) polyp respiration.
Colony respiration:
colony.R.effects <- predictorEffect(c("Org.Number"), Rmodel2) #set x axis
colony.effR<-plot(colony.R.effects,
lines=list(multiline=TRUE, #color lines
col=c("purple"),
lty=c(1)),
confint=list(style="bands", alpha=0.3), #set conf int
lwd=4,
axes=list(y=list(lim=c(-.4, 0.1))),
type="response", #set response scale
ylab=expression(bold("Colony Total Respiration")),
legend.position="top",
xlab=expression(bold("Number of Fused Polyps")),
main="");colony.effR
Polyp specific respiration:
polyp.R.effects <- predictorEffect(c("Org.Number"), Rmodel3) #set x axis
polyp.effR<-plot(polyp.R.effects,
lines=list(multiline=TRUE, #color lines
col=c("orange"),
lty=c(1)),
confint=list(style="bands", alpha=0.3), #set conf int
lwd=4,
axes=list(y=list(lim=c(-.4, 0.1))),
type="response", #set response scale
ylab=expression(bold("Polyp Specific Respiration")),
legend.position="top",
xlab=expression(bold("Number of Fused Polyps")),
main="");polyp.effR
#lattice=list(key.args=list(space="right",
#border=FALSE,
#title=expression(bold("Temperature - Fusion")),
#cex=1,
#cex.title=1)))
In fused colonies, polyp-specific respiration and total colony respiration are both significantly affected by the size of the colony (number of fused polyps). Polyp specific respiration seems to decrease in larger colonies while total colony respiration increases. This may suggest that each constituent polyp “saves” energy being part of a larger fused colony. This is very interesting to think about!
Analyze net photosynthesis at the level of the whole colony, or normalized at polyp-specific rates within each colony. Analyze across total number of polyps fused together.
Build linear mixed effect model and examine for Net Photosynthesis at the colony level.
Pmodel2<-lmer((P.nmol.org.min*Org.Number)~Org.Number + (1|Run/dpf), data=fused_oxygen) #run nested within day
anova(Pmodel2, type="II")
## Type II Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Org.Number 0.30307 0.30307 1 18 30.837 2.851e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Check assumptions of model for residual normality and variance. Assumptions are met.
qqPlot(residuals(Pmodel2))
## 38 31
## 18 11
hist(residuals(Pmodel2))
leveneTest(residuals(Pmodel2)~as.factor(Org.Number), data=fused_oxygen)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 7 1.1103 0.4161
## 12
Build linear mixed effect model and examine for Net Photosynthesis at the polyp level.
Pmodel3<-lmer(P.nmol.org.min~Org.Number + (1|Run/dpf), data=fused_oxygen) #run nested within day
anova(Pmodel3, type="II")
## Type II Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Org.Number 0.0054339 0.0054339 1 18 5.1573 0.03566 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Check assumptions of model for residual normality and variance. Assumptions met.
qqPlot(residuals(Pmodel3))
## 27 21
## 7 1
hist(residuals(Pmodel3))
leveneTest(residuals(Pmodel3)~as.factor(Org.Number), data=fused_oxygen)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 7 0.5616 0.7734
## 12
Display linear regression of 1) colony net photosynthesis and 2) polyp net photosynthesis.
Colony net photosynthesis:
colony.NP.effects <- predictorEffect(c("Org.Number"), Pmodel2) #set x axis
colony.effNP<-plot(colony.NP.effects,
lines=list(multiline=TRUE, #color lines
col=c("purple"),
lty=c(1)),
confint=list(style="bands", alpha=0.3), #set conf int
lwd=4,
axes=list(y=list(lim=c(-0.1, 0.8))),
type="response", #set response scale
ylab=expression(bold("Colony Total Net Photosynthesis")),
legend.position="top",
xlab=expression(bold("Number of Fused Polyps")),
main="");colony.effNP
Polyp net photosynthesis:
polyp.NP.effects <- predictorEffect(c("Org.Number"), Pmodel3) #set x axis
polyp.effNP<-plot(polyp.NP.effects,
lines=list(multiline=TRUE, #color lines
col=c("orange"),
lty=c(1)),
confint=list(style="bands", alpha=0.3), #set conf int
lwd=4,
axes=list(y=list(lim=c(-0.1, 0.8))),
type="response", #set response scale
ylab=expression(bold("Polyp Specific Net Photosynthesis")),
legend.position="top",
xlab=expression(bold("Number of Fused Polyps")),
main="");polyp.effNP
In fused colonies, net photosynthesis increases with larger fused colonies (significant), but polyp-specific net photosynthesis does not increase (not significant). This suggests that the capacity of each polyp to photosynthesize is not enhanced with fusion, but rather there is an additive effect of photosynthesis across the larger colony.
Analyze gross photosynthesis at the level of the whole colony, or normalized at polyp-specific rates within each colony. Analyze across total number of polyps fused together.
Build linear mixed effect model and examine for Gross Photosynthesis at the colony level.
GPmodel2<-lmer((GP.nmol.org.min*Org.Number)~Org.Number + (1|Run/dpf), data=fused_oxygen) #run nested within day
anova(GPmodel2, type="II")
## Type II Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Org.Number 0.55445 0.55445 1 18 28.939 4.114e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Check assumptions of model for residual normality and variance. Assumptions are met.
qqPlot(residuals(GPmodel2))
## 38 36
## 18 16
hist(residuals(GPmodel2))
leveneTest(residuals(GPmodel2)~as.factor(Org.Number), data=fused_oxygen)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 7 1.5169 0.2509
## 12
Build linear mixed effect model and examine for Gross Photosynthesis at the polyp level.
GPmodel3<-lmer(GP.nmol.org.min~Org.Number + (1|Run/dpf), data=fused_oxygen) #run nested within day
anova(GPmodel3, type="II")
## Type II Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Org.Number 0.013236 0.013236 1 18 8.0413 0.01096 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Check assumptions of model for residual normality and variance. Potential normality violation, revisit.
qqPlot(residuals(GPmodel3))
## 38 29
## 18 9
hist(residuals(GPmodel3))
leveneTest(residuals(GPmodel3)~as.factor(Org.Number), data=fused_oxygen)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 7 0.9294 0.5182
## 12
Display linear regression of 1) colony respiration and 2) polyp respiration.
Colony gross photosynthesis:
colony.GP.effects <- predictorEffect(c("Org.Number"), GPmodel2) #set x axis
colony.effGP<-plot(colony.GP.effects,
lines=list(multiline=TRUE, #color lines
col=c("purple"),
lty=c(1)),
confint=list(style="bands", alpha=0.3), #set conf int
lwd=4,
axes=list(y=list(lim=c(-0.1, 0.8))),
type="response", #set response scale
ylab=expression(bold("Colony Total Gross Photosynthesis")),
legend.position="top",
xlab=expression(bold("Number of Fused Polyps")),
main="");colony.effGP
Polyp gross photosynthesis:
polyp.GP.effects <- predictorEffect(c("Org.Number"), GPmodel3) #set x axis
polyp.effGP<-plot(polyp.GP.effects,
lines=list(multiline=TRUE, #color lines
col=c("orange"),
lty=c(1)),
confint=list(style="bands", alpha=0.3), #set conf int
lwd=4,
axes=list(y=list(lim=c(-0.1, 0.8))),
type="response", #set response scale
ylab=expression(bold("Polyp Specific Gross Photosynthesis")),
legend.position="top",
xlab=expression(bold("Number of Fused Polyps")),
main="");polyp.effGP
In fused colonies, gross photosynthesis has similar results as net photosynthesis, which makes sense. Gross photosynthesis increases with larger fused colonies (significant), but polyp-specific net photosynthesis does not increase (not significant). This suggests that the capacity of each polyp to photosynthesize is not enhanced with fusion, but rather there is an additive effect of photosynthesis across the larger colony.
Build linear mixed effect model and examine for P:R ratio (gross photosynthesis : respiration). Ratios are the same when calculated for colony or polyp level, so displaying and analyzing at the polyp level here.
PRmodel2<-lmer(ratio~Org.Number + (1|Run/dpf), data=fused_oxygen)
anova(PRmodel2, type="II")
## Type II Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Org.Number 1.0467 1.0467 1 18 0.9065 0.3537
Check assumptions of model for residual normality and variance. Assumptions met.
qqPlot(residuals(PRmodel2))
## 39 22
## 19 2
hist(residuals(PRmodel2))
leveneTest(residuals(PRmodel2)~as.factor(Org.Number), data=fused_oxygen)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 7 3.9403 0.01833 *
## 12
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
P:R ratios:
ratio.effects <- predictorEffect(c("Org.Number"), PRmodel2) #set x axis
ratio.effGP<-plot(ratio.effects,
lines=list(multiline=TRUE, #color lines
col=c("black"),
lty=c(1)),
confint=list(style="bands", alpha=0.3), #set conf int
lwd=4,
axes=list(y=list(lim=c(-0.1, 8))),
type="response", #set response scale
ylab=expression(bold("P:R Ratio")),
legend.position="top",
xlab=expression(bold("Number of Fused Polyps")),
main="");ratio.effGP
P:R ratio does not significantly change over colony size. Juveniles at this age maintain positive P:R ratios (above 1) regardless of size. This is likely achieved by reduced polyp-specific respiration balanced by increased colony total respiration, keeping respiratory demand consistent across colony sizes.
Combine regression plts above into a single panel organized by metric.
fused_regressions<-plot_grid(polyp.effR, colony.effR, polyp.effNP, colony.effNP, polyp.effGP, colony.effGP, ratio.effGP, labels = c("p=0.01", "p<0.01", "p>0.05", "p=0.04", "p>0.05", "p=0.03", "p>0.05"), label_size = 16, ncol=2, nrow=4, rel_heights= c(1), rel_widths = c(1), align="h")
ggsave(filename="Mcap2020/Figures/Respiration/fused_regressions.pdf", plot=fused_regressions, dpi=300, width=10, height=16, units="in")
In this data, we saw a few things that are particularly interesting to think about:
(1) Respiration and photosynthesis are different between lifestages. Larval respiration and photosynthesis are low in larvae, and just maintain energetic balance (P meeting R demands).
(2) Respiration increases at recruitment (2-4 days after recruitment, 8-10 days after fertilization) and is not met by photosynthesis, highlighting energetic vulnerabilities at this stage.
(3) For older juveniles (40 days post fertilization), there is excess in photosynthesis to meet respiratory demand that is available for growth and calcification.
(4) In juveniles, those that are fused have higher respiration and photosynthesis as the size of the colony (number of polyps fused together) is larger. Photosynthesis meets respiratory demand and is in excess, providing energy for growth and calcification at this stage.
(5) Interestingly, polyp-specific respiration decreases in larger fused colonies while total respiration increases. This suggests that by joining together in a larger colony, each constituent polyp may “save” energy.
(6) Raises interesting questions for integration timing of symbiosis after settlement.